//============================================================================
// Name        : FabVoronoi.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <iterator>
#include <iomanip>
#include "Topol.h"
#include "FAtoms.h"
#include "HeaderTrj.h"
#include "ResidueCM.h"
#include "FVoronoi.h"

using namespace Topol_NS;
using namespace std;
int main(int argc, char ** argv)
{
	bool bHyd=true;
	int narg=1;
	string filein,fileout,filepdb;
	ifstream fin,fpdb;
	ofstream fout;
	string comm0(argv[0]);
	size_t mypos=comm0.find_last_of("/")+1;
	size_t length=comm0.size()-mypos;
	string command=comm0.substr(mypos,length);
	const string Usage=string("Usage:\t"+ command + "-b startframe[0] -e endframe[end] " +
			"-skip skipfrequency[1] \n\t -i filein -o fileout -pdb filepdb -[no]hyd \n "
			"\n\t Default values in square backets []\n");
	try{
		if(argc == 1) throw Usage;
	}
	catch(const string & s){
		cout << s << endl;
		exit(1);
	}

	int nstart=0;
	int nend=-1;
	int nskip=1;
// input read
	do{
		string myarg0;

		myarg0.assign(argv[narg]);
		try{
			if(myarg0.find("-i") != string::npos){
				narg++;
				if(narg == argc ||argv[narg][0] == '-') throw string(" filename expected ");
				filein.assign(argv[narg]);
				fin.open(filein.c_str(),ios::in);
			}
			else if(myarg0.find("-o") != string::npos){
				narg++;

				if(narg == argc || argv[narg][0] == '-') throw string(" filename expected ");
				fileout.assign(argv[narg]);
				fout.open(fileout.c_str(),ios::out);
			}
			else if(myarg0.find("-b") != string::npos){
				narg++;
				if(narg == argc || argv[narg][0] == '-') throw string("begin time ");
				stringstream(argv[narg])>> nstart;
			}
			else if(myarg0.find("-e") != string::npos){
				narg++;
				if(narg == argc || argv[narg][0] == '-') throw string("end time ");
				stringstream(argv[narg])>> nend;
			}
			else if(myarg0.find("-skip") != string::npos){
				narg++;
				if(narg == argc || argv[narg][0] == '-') throw string("skip time frequency");
				stringstream(argv[narg])>> nskip;
			}
			else if(myarg0.find("-pdb") != string::npos){
				narg++;
				if(narg == argc ||argv[narg][0] == '-') throw string(" filename expected ");
				filepdb.assign(argv[narg]);
				fpdb.open(filepdb.c_str(),ios::in);
			}
			else if(myarg0.find("-nohyd") != string::npos){
				bHyd=false;
			}
			else if(myarg0.find("-hyd") != string::npos){
				bHyd=true;
			}
			else{
				throw Usage;
			}
		}
		catch(const string  & s){
			cout << s << endl;
			exit(1);
		};

	} while(++narg<argc);

// read pdb file to construct topology
	vector<string> data;
	for(string str;getline(fpdb,str);){
		data.push_back(str);
	}
	Topol MyTop(data);

	int natoms=MyTop.Size();
	int nresid=MyTop.ResSize();
	FAtoms atm;
	atm.setDim(natoms);

// Map atoms to residues to use in Voronoi. Eliminate hydrogens if needed
	vector<int> * resmap=new vector<int> [nresid];
    string * ResLabels=new string [nresid];
    string label;

    for(int i=0;i<nresid;i++)
    	ResLabels[i]=MyTop.ResInfo(i);

    for(int i=0;i<natoms;i++){
    	int j=MyTop.ResIndex(i);
    	if(bHyd) resmap[j].push_back(i);
    	else {
    		label=MyTop.AtomInfo(i);
    		if(label.compare(0,1,"H",0,1) != 0) resmap[j].push_back(i);
    	}
    }

    ResidueCM MyCM;
    MyCM(natoms,nresid,resmap,ResLabels);
    AtomIndex cidx[TYPES];
	for(int i=0;i<TYPES;i++){
		vector<int> ty=MyTop.getCidx(i);
		cidx[i].set(ty.size(),&ty[0]);
	}

	MyCM.ResidueSelect(TYPES,cidx);
// Setup Voronoi

	FVoronoi * vor=new FVoronoi(natoms,MyTop.getAtomName(),bHyd);
	vor->setTypes(TYPES,cidx,MyTop.getTypeNames());

// Read header of dcd file
	fin.seekg(0,ios::end);
	int len=fin.tellg();
	fin.seekg(0,ios::beg);

	HeaderTrj oppo;
	fin>>oppo;
	try{
		if(!oppo.check(MyTop.Size())) throw string("Number of atoms in the pdb and trajectory files does not match!");
	}
	catch(const string & s){
		cout << s<<endl;
		exit(-1);
	}
// Read trajectory and compute Voronoi staff
	int ntime=0;
	while(fin.tellg() < len){
		fin>> atm;
		if(ntime < nstart) {ntime++;continue;}
		if(nend != -1 && ntime > nend) break;
		if((ntime-nstart)%nskip) {ntime++;continue;}
		cout << fixed << setw(5) << "----> Step No. " << ntime <<"\n";
		vor->Start(static_cast<float> (ntime),atm);
		vor->getData();
// to be done:  write volumes according to Fabio's fortran code
		fout << *vor;
		ntime++;
	};
	cout << "\nProgram completed: Output data written to " + fileout << "\n\n";
/*
	vector<string> data((istream_iterator<string>(fin)),istream_iterator<string>());
	for(int i=0;i<data.size();i++) cout << data[i] << endl;
	*/
    return 0;
}   //end of main function
